# Title: Public Trust and Collaborative Governance: An Instrumental Variable Approach
# Yixin Liu 2022
# Public Management Review
# PMR Replication file for Study 2

library(ggplot2) # plot
library(ggridges) # overlapping line plots
library(pwr) # power
library(ggpubr) # plot
library(effsize) # effect size
library(multiwayvcov) # cluster se
library(lmtest) # cluster se
library(egg) # ggplot(theme_article)
library(cowplot) # ggplot output
library(Rmisc) # plot stat table function
library(stargazer) # latex
library(xtable) # latex
library(texreg) # latex
library(ivreg)
library(standardize)
library(Rmisc) # summaryse function

####Function####
element_textbox <- function(...) {
  el <- element_text(...)
  class(el) <- c("element_textbox", class(el))
  el
}

element_grob.element_textbox <- function(element, ...) {
  text_grob <- NextMethod()
  rect_grob <- element_grob(calc_element("strip.background", theme_bw()))
  
  ggplot2:::absoluteGrob(
    grid::gList(
      element_grob(calc_element("strip.background", theme_bw())),
      text_grob
    ),
    height = grid::grobHeight(text_grob), 
    width = grid::unit(1, "npc")
  )
}

# Setup working directory
WD = getwd()
setwd(WD)

####Data####
exp<-read.csv("study2_data.csv")

####Trust####
# combine measurement
exp$trust_combine<-(exp$competence+exp$benevolence+exp$honesty)/3

exp$trustprime<-NA
exp$trustprime[exp$rand=="control"]<-1
exp$trustprime[exp$rand=="lt"]<-0
exp$trustprime[exp$rand=="ht"]<-2

exp$rand[exp$rand=="control"]<-"Control"
exp$rand[exp$rand=="lt"]<-"Corrupt"
exp$rand[exp$rand=="ht"]<-"Honest"

trust.aov<-aov(trust_combine~rand, data = exp)
summary(trust.aov)

####Figure 3####
# left panel
trust.df<-summarySE(exp, measurevar="trust_combine", 
                    groupvars=c("rand"),
                    conf.interval = 0.95,
                    na.rm = T)

trust.df$rand<-as.factor(trust.df$rand)
trust.df$rand<-ordered(trust.df$rand,level=c("Corrupt","Control","Honest"))

trust.plot<-ggplot(trust.df,aes(x=rand, y=trust_combine, group=1)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(size=3)) +
  geom_pointrange(aes(ymin=trust_combine-ci, ymax=trust_combine+ci)) +
  #scale_y_continuous(lim=c(30,80)) +
  labs(x = element_blank(),
       y = element_blank(),
       title = "Trust in Government")+
  theme_article()

trust.plot<-trust.plot + theme(legend.position = "none",
                               axis.title.y = element_text(size = 12),
                               axis.text.x=element_text(size = 12),
                               axis.text.y=element_text(size = 10),
                               strip.background = element_blank(),
                               strip.text.x = element_blank(),
                               axis.ticks.x = element_blank(),
                               plot.title = element_textbox(size = 16,
                                                            hjust = 0.5, margin = margin(t = 5, b = 5)))

trust.plot

# right panel
private.df<-summarySE(exp, measurevar="PPP", 
                      groupvars=c("rand"),
                      conf.interval = 0.95,
                      na.rm = T)

community.df<-summarySE(exp, measurevar="PCP", 
                        groupvars=c("rand"),
                        conf.interval = 0.95,
                        na.rm = T)

coproduce.df<-summarySE(exp, measurevar="coproduce", 
                        groupvars=c("rand"),
                        conf.interval = 0.95,
                        na.rm = T)

names(private.df)[names(private.df) == "PPP"] <- "dv"
private.df$group<-"Public-Private Partnerships"

names(community.df)[names(community.df) == "PCP"] <- "dv"
community.df$group<-"Public-Citizen Partnerships"

names(coproduce.df)[names(coproduce.df) == "coproduce"] <- "dv"
coproduce.df$group<-"Willingness to Coproduce"

df.dv<-rbind(private.df,community.df,coproduce.df)
df.dv$rand<-ordered(df.dv$rand,level=c("Corrupt","Control","Honest"))

df.dv$group <- factor(df.dv$group, 
                      levels = c("Public-Citizen Partnerships", 
                                 "Willingness to Coproduce",
                                 "Public-Private Partnerships"))

dv.plot<-ggplot(df.dv,aes(x=rand, y=dv, group=group, shape=group)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(shape=group), size = 4) +
  geom_pointrange(aes(ymin=dv-ci, ymax=dv+ci)) +
  #scale_colour_manual(values = c("#00BFC4", "#F8766D", "#7CAE00")) +
  scale_y_continuous(lim=c(50,85)) +
  labs(x = element_blank(),
       y = element_blank(),
       title = "Percieved Collaborative Governance")+
  theme_article()

dv.plot<-dv.plot + theme(legend.position = c(0.75, 0.15),
                         legend.direction = "vertical",
                         legend.title = element_blank(),
                         legend.background = element_blank(),
                         legend.text = element_text(size = 12),
                         strip.background = element_blank(),
                         strip.text.x = element_blank(),
                         axis.ticks.x = element_blank(),
                         axis.title.y = element_text(size = 12),
                         axis.text.x=element_text(size = 12),
                         axis.text.y=element_text(size = 10),
                         plot.title = element_textbox(size = 16,
                                                      hjust = 0.5, margin = margin(t = 5, b = 5)))

dv.plot

####trust + dv####
combine.plot<-ggarrange(trust.plot,
                        dv.plot,
                        ncol = 2, nrow = 1)
combine.plot<-plot_grid(combine.plot)
combine.plot
#ggsave(file = "figure3.jpg", width = 12, height = 5, dpi = 666)

####Table 5####
# stage 1 (model 1)
trust1<-lm(trust_combine~trustprime,data=exp)
summary(trust1)

# stage 2 (model 2)
iv1 = ivreg(PPP~trust_combine
            | trustprime, data = exp)
summary(iv1)

# stage 2 (model 3)
iv2 = ivreg(PCP~trust_combine
            | trustprime, data = exp)
summary(iv2)

# stage 2 (model 4)
iv3 = ivreg(coproduce~trust_combine 
            | trustprime, data = exp)
summary(iv3)

# regression table
main1<-list(trust1,iv1,iv2,iv3)
texreg(main1,
       custom.model.names=c("(1)","(2)","(3)","(4)"),
       custom.coef.map=list("trust_combine"="Trust in Government",
                            "trustprime"="Government Integrity",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. Government integrity is measured 
                           by coding experimental conditions as: 
                           1 = Honest; 0.5 = Control; 0 = Corrupt."),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       caption.above=TRUE,digits=2,dcolumn=TRUE,
       caption="2SLS Estimate in Study 2")

####Table 6####
# subgroup separate instrument
exp.honest<-subset(exp, exp$trustprime!=0)
exp.honest$trustprime<-ifelse(exp.honest$trustprime==1,0,1)
exp.corrupt<-subset(exp, exp$trustprime!=2)
exp.corrupt$trustprime<-ifelse(exp.corrupt$trustprime==1,0,1)

# Panel A
# stage 1 (model A1)
honest1<-lm(trust_combine~trustprime,data=exp.honest)
summary(honest1)

# stage 2 (model A2)
iv1.honest = ivreg(PPP~trust_combine
                   | trustprime, data = exp.honest)
summary(iv1.honest)

# stage 2 (model A3)
iv2.honest = ivreg(PCP~trust_combine
                   | trustprime, data = exp.honest)
summary(iv2.honest)

# stage 2 (model A4)
iv3.honest = ivreg(coproduce~trust_combine 
                   | trustprime, data = exp.honest)
summary(iv3.honest)

# Panel B
# stage 1 (model B1)
corrupt1<-lm(trust_combine~trustprime,data=exp.corrupt)
summary(corrupt1)

# stage 2 (model B2)
iv1.corrupt = ivreg(PPP~trust_combine
                    | trustprime, data = exp.corrupt)
summary(iv1.corrupt)

# stage 2 (model B3)
iv2.corrupt = ivreg(PCP~trust_combine
                    | trustprime, data = exp.corrupt)
summary(iv2.corrupt)

# stage 2 (model B4)
iv3.corrupt = ivreg(coproduce~trust_combine 
                    | trustprime, data = exp.corrupt)
summary(iv3.corrupt)

# regression table
# Panel A
honest1<-list(honest1,iv1.honest,iv2.honest,iv3.honest)
texreg(honest1,
       custom.model.names=c("(1)","(2)","(3)","(4)"),
       custom.coef.map=list("trust_combine"="Trust in Government",
                            "trustprime"="Honest",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. Honest is measured by 
                            coding experimental conditions as: 1 = Honest; 0 = Control."),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       caption.above=TRUE,digits=2,dcolumn=TRUE,
       caption="2SLS Estimate Using Honest and Corrupt as Separate Instruments in Study 2 (Panel A)")

# Panel B
corrupt1<-list(corrupt1,iv1.corrupt,iv2.corrupt,iv3.corrupt)
texreg(corrupt1,
       custom.model.names=c("(1)","(2)","(3)","(4)"),
       custom.coef.map=list("trust_combine"="Trust in Government",
                            "trustprime"="Corrupt",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. Corrupt is measured by 
                            coding experimental conditions as: 1 = Corrupt; 0 = Control."),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       caption.above=TRUE,digits=2,dcolumn=TRUE,
       caption="2SLS Estimate Using Honest and Corrupt as Separate Instruments in Study 2 (Panel B)")

####Appendix B####
#####Appendix B.2####
# study 2 sample
vars<-data.frame(exp$sex,
                 exp$white,exp$black,exp$hispanic,exp$asian,exp$other,
                 exp$ideology,
                 exp$age,
                 exp$education,
                 exp$income,
                 exp$gw,exp$gwh,exp$climatebelief,
                 exp$covid_concern,exp$covid_job)

Mean=sapply(vars, mean, na.rm = TRUE)
SD=sapply(vars, sd, na.rm = TRUE)
Min=sapply(vars, min, na.rm = TRUE)
Max=sapply(vars, max, na.rm = TRUE)

female.aov <- aov(sex ~ rand, data = exp)
summary(female.aov)
female<-summary(female.aov)[[1]][,5][1]

wht.aov <- aov(white ~ rand, data = exp)
summary(wht.aov)
wht<-summary(wht.aov)[[1]][,5][1]
blk.aov <- aov(black ~ rand, data = exp)
summary(blk.aov)
blk<-summary(blk.aov)[[1]][,5][1]
his.aov <- aov(hispanic ~ rand, data = exp)
summary(his.aov)
his<-summary(his.aov)[[1]][,5][1]
asian.aov <- aov(asian ~ rand, data = exp)
summary(asian.aov)
asian<-summary(asian.aov)[[1]][,5][1]
other.aov <- aov(other ~ rand, data = exp)
summary(other.aov)
other<-summary(other.aov)[[1]][,5][1]

ideo.aov <- aov(ideology ~ rand, data = exp)
summary(ideo.aov)
ideo<-summary(ideo.aov)[[1]][,5][1]

age.aov <- aov(age ~ rand, data = exp)
summary(age.aov)
age<-summary(age.aov)[[1]][,5][1]

educ.aov <- aov(education ~ rand, data = exp)
summary(educ.aov)
educ<-summary(educ.aov)[[1]][,5][1]

inc.aov <- aov(income ~ rand, data = exp)
summary(inc.aov)
inc<-summary(inc.aov)[[1]][,5][1]

gw.aov <- aov(gw ~ rand, data = exp)
summary(gw.aov)
gw<-summary(gw.aov)[[1]][,5][1]

gwh.aov <- aov(gwh ~ rand, data = exp)
summary(gwh.aov)
gwh<-summary(gwh.aov)[[1]][,5][1]

climate.aov <- aov(climatebelief ~ rand, data = exp)
summary(climate.aov)
climate<-summary(climate.aov)[[1]][,5][1]

covid.aov <- aov(covid_concern ~ rand, data = exp)
summary(covid.aov)
covid<-summary(covid.aov)[[1]][,5][1]

covid.job.aov <- aov(covid_job ~ rand, data = exp)
summary(covid.job.aov)
covidjob<-summary(covid.job.aov)[[1]][,5][1]

# sample table
des_stat<-matrix(NA,length(Mean),5)
colnames(des_stat)<-c("Mean","SD","Min","Max","Randomization Check (P-value)")
rownames(des_stat)<-c("Female","White","Black","Hispanic","Asian","Other",
                      "Ideology","Age","Education","Income", "GW", 
                      "GWH", "Climate Change Beliefs",
                      "COVID Concern", "COVID Job Concern")

des_stat[,1]<-Mean
des_stat[,2]<-SD
des_stat[,3]<-Min
des_stat[,4]<-Max
des_stat[,5]<-c(female,wht,blk,his,asian,other,ideo,age,
                educ,inc,gw,gwh,climate,covid,covidjob)

des_stat<-round(des_stat,digits = 2)

# write table
xtable(des_stat, digits=2)

#####Appendix B.3####
# Perceived Competence, Benevolence, and Honest by Experimental Conditions
# competence
competence.df<-summarySE(exp, measurevar="competence", 
                         groupvars=c("rand"),
                         conf.interval = 0.95,
                         na.rm = T)

competence.df$rand<-as.factor(competence.df$rand)
competence.df$rand<-ordered(competence.df$rand,level=c("Corrupt","Control","Honest"))

competence.plot<-ggplot(competence.df,aes(x=rand, y=competence, group=1)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(size=3)) +
  geom_pointrange(aes(ymin=competence-ci, ymax=competence+ci)) +
  scale_y_continuous(lim=c(3, 5.5)) +
  labs(x = element_blank(),
       y = element_blank(),
       title = "Competence")+
  theme_article()

competence.plot<-competence.plot + theme(legend.position = "none",
                                         axis.title.y = element_text(size = 12),
                                         axis.text.x=element_text(size = 12),
                                         axis.text.y=element_text(size = 10),
                                         strip.background = element_blank(),
                                         strip.text.x = element_blank(),
                                         axis.ticks.x = element_blank(),
                                         plot.title = element_textbox(
                                           hjust = 0.5, margin = margin(t = 5, b = 5)))

competence.plot

# benevolence
benevolence.df<-summarySE(exp, measurevar="benevolence", 
                          groupvars=c("rand"),
                          conf.interval = 0.95,
                          na.rm = T)

benevolence.df$rand<-as.factor(benevolence.df$rand)
benevolence.df$rand<-ordered(benevolence.df$rand,level=c("Corrupt","Control","Honest"))

benevolence.plot<-ggplot(benevolence.df,aes(x=rand, y=benevolence, group=1)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(size=3)) +
  geom_pointrange(aes(ymin=benevolence-ci, ymax=benevolence+ci)) +
  scale_y_continuous(lim=c(3,5.5)) +
  labs(x = element_blank(),
       y = element_blank(),
       title = "Benevolence")+
  theme_article()

benevolence.plot<-benevolence.plot + theme(legend.position = "none",
                                           axis.title.y = element_text(size = 12),
                                           axis.text.x=element_text(size = 12),
                                           axis.text.y=element_text(size = 10),
                                           strip.background = element_blank(),
                                           strip.text.x = element_blank(),
                                           axis.ticks.x = element_blank(),
                                           plot.title = element_textbox(
                                             hjust = 0.5, margin = margin(t = 5, b = 5)))

benevolence.plot

# honesty
honesty.df<-summarySE(exp, measurevar="honesty", 
                      groupvars=c("rand"),
                      conf.interval = 0.95,
                      na.rm = T)

honesty.df$rand<-as.factor(honesty.df$rand)
honesty.df$rand<-ordered(honesty.df$rand,level=c("Corrupt","Control","Honest"))

honesty.plot<-ggplot(honesty.df,aes(x=rand, y=honesty, group=1)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(size=3)) +
  geom_pointrange(aes(ymin=honesty-ci, ymax=honesty+ci)) +
  scale_y_continuous(lim=c(3,5.5)) +
  labs(x = element_blank(),
       y = element_blank(),
       title = "Honesty")+
  theme_article()

honesty.plot<-honesty.plot + theme(legend.position = "none",
                                   axis.title.y = element_text(size = 12),
                                   axis.text.x=element_text(size = 12),
                                   axis.text.y=element_text(size = 10),
                                   strip.background = element_blank(),
                                   strip.text.x = element_blank(),
                                   axis.ticks.x = element_blank(),
                                   plot.title = element_textbox(
                                     hjust = 0.5, margin = margin(t = 5, b = 5)))

honesty.plot

# plot combine
trust_combine.plot<-ggarrange(competence.plot,
                              benevolence.plot+ theme(axis.text.y=element_blank(),
                                                      axis.ticks.y = element_blank()),
                              honesty.plot+ theme(axis.text.y=element_blank(),
                                                  axis.ticks.y = element_blank()),
                              ncol = 3, nrow = 1)
trust_combine.plot<-plot_grid(trust_combine.plot)
trust_combine.plot
#ggsave(file = "b3.jpg", width = 8, height = 5, dpi = 666)

# ANOVA
competence.aov<-aov(competence~rand, data = exp)
summary(competence.aov)

benevolence.aov<-aov(benevolence~rand, data = exp)
summary(benevolence.aov)

honesty.aov<-aov(honesty~rand, data = exp)
summary(honesty.aov)

#####Appendix B.4####
trustcat<-lm(trust_combine~rand,data=exp)
summary(trustcat)

cat1<-lm(PPP~rand, data = exp)
summary(cat1)

cat2<-lm(PCP~rand, data = exp)
summary(cat2)

cat3<-lm(coproduce~rand, data = exp)
summary(cat3)

# regression table
cat<-list(trustcat,cat1,cat2,cat3)
texreg(cat,
       custom.model.names=c("(1)","(2)","(3)","(4)"),
       custom.coef.map=list("randHonest"="Honest",
                            "randCorrupt"="Corrupt",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. "),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       caption.above=TRUE,digits=2,dcolumn=TRUE)

#####Appendix B.5####
covid1 = lm(trust_combine~covid_concern+covid_job, data = exp)
summary(covid1)

covid2 = ivreg(coproduce~trust_combine*covid_concern
               | trustprime*covid_concern, data = exp)
summary(covid2)

covid3 = ivreg(coproduce~trust_combine*covid_job
               | trustprime*covid_job, data = exp)
summary(covid3)

# regression model
covid<-list(covid1,covid2,covid3)
texreg(covid,
       custom.model.names=c("(1)","(2)","(3)"),
       custom.coef.map=list("trust_combine"="Trust in Government",
                            "covid_concern"="Covid Concern",
                            "covid_job"="Covid Job Concern",
                            "trust_combine:covid_concern"="Trust in Government$\\times$Covid Concern",
                            "trust_combine:covid_job"="Trust in Government$\\times$Covid Job Concern",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. "),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       caption.above=TRUE,digits=2,dcolumn=TRUE)

